{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "from IPython.display import display, Image\n",
    "\n",
    "from rmgpy.tools.uncertainty import Uncertainty, process_local_results\n",
    "from rmgpy.tools.canteraModel import get_rmg_species_from_user_species\n",
    "from rmgpy.species import Species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First Order Local Uncertainty Analysis for Chemical Reaction Systems\n",
    "\n",
    "This IPython notebook performs first order local uncertainty analysis for a chemical reaction system\n",
    "using a RMG-generated model.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Define mechanism files and simulation settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a small ethane pyrolysis model\n",
    "\n",
    "# Must use annotated chemkin file\n",
    "chemkin_file = 'data/parse_source/chem_annotated.inp'\n",
    "dict_file = 'data/parse_source/species_dictionary.txt'\n",
    "\n",
    "# Initialize the Uncertainty class instance and load the model\n",
    "uncertainty = Uncertainty(output_directory='./temp/uncertainty')\n",
    "uncertainty.load_model(chemkin_file, dict_file)\n",
    "\n",
    "# Map the species to the objects within the Uncertainty class\n",
    "ethane = Species().from_smiles('CC')\n",
    "C2H4 = Species().from_smiles('C=C')\n",
    "mapping = get_rmg_species_from_user_species([ethane, C2H4], uncertainty.species_list)\n",
    "\n",
    "# Define the reaction conditions\n",
    "initial_mole_fractions = {mapping[ethane]: 1.0}\n",
    "T = (1300, 'K')\n",
    "P = (1, 'bar')\n",
    "termination_time = (0.5, 'ms')\n",
    "sensitive_species=[mapping[ethane], mapping[C2H4]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Run sensitivity analysis\n",
    "\n",
    "Local uncertainty analysis uses the results from a first-order sensitivity analysis. This analysis is done using RMG's native solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform the sensitivity analysis\n",
    "uncertainty.sensitivity_analysis(initial_mole_fractions, sensitive_species, T, P, termination_time, number=5, fileformat='.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Show the sensitivity plots\n",
    "for species in sensitive_species:\n",
    "    print('{}: Reaction Sensitivities'.format(species))\n",
    "    index = species.index\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory,'solver','sensitivity_1_SPC_{}_reactions.png'.format(index))))\n",
    "    \n",
    "    print('{}: Thermo Sensitivities'.format(species))\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory,'solver','sensitivity_1_SPC_{}_thermo.png'.format(index))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Uncertainty assignment and propagation of uncorrelated parameters\n",
    "\n",
    "If we want to run local uncertainty analysis, we must assign all the uncertainties using the `Uncertainty` class' `assignParameterUncertainties` function. `ThermoParameterUncertainty` and `KineticParameterUncertainty` classes may be customized and passed into this function if non-default constants for constructing the uncertainties are desired. This must be done after the parameter sources are properly extracted from the model.\n",
    "\n",
    "### Thermo Uncertainty\n",
    "\n",
    "Each species is assigned a uniform uncertainty distribution in free energy:\n",
    "\n",
    "$$G \\in [G_{min},G_{max}]$$\n",
    "\n",
    "We will propogate the standard deviation in free energy, which for a uniform distribution is defined as follows:\n",
    "\n",
    "$$\\Delta G = \\frac{1}{\\sqrt{12}}(G_{max} - G_{min})$$\n",
    "\n",
    "Several parameters are used to formulate $\\Delta G$.  These are $\\Delta G_\\mathrm{library}$, $\\Delta G_\\mathrm{QM}$, $\\Delta G_\\mathrm{GAV}$, and $\\Delta _\\mathrm{group}$.\n",
    "        \n",
    "$$\\Delta G = \\delta_\\mathrm{library} \\Delta G_\\mathrm{library} + \\delta_\\mathrm{QM} \\Delta G_\\mathrm{QM} + \\delta_\\mathrm{GAV} \\left( \\Delta G_\\mathrm{GAV} + \\sum_{\\mathrm{group}\\; j} d_{j} \\Delta G_{\\mathrm{group},j} \\right)$$\n",
    "\n",
    "where $\\delta$ is the Kronecker delta function which equals one if the species thermochemistry parameter contains the particular source type and $d_{j}$ is the degeneracy (number of appearances) of the thermo group used to construct the species thermochemistry in the group additivity method.\n",
    "\n",
    "### Kinetics Uncertainty\n",
    "\n",
    "Each reaction is assigned a uniform uncertainty distribution in the overall $\\ln k$, or $\\ln A$:\n",
    "\n",
    "$$\\ln k \\in [\\ln(k_{min}),\\ln(k_{max})]$$\n",
    "\n",
    "Again, we use the standard deviation of this distribution:\n",
    "\n",
    "$$\\Delta \\ln(k) = \\frac{1}{\\sqrt{12}}(\\ln k_{max} - \\ln k_{min})$$\n",
    "\n",
    "The parameters used to formulate $\\Delta  \\ln k$ are $\\Delta \\ln k_\\mathrm{library}$, $\\Delta \\ln k_\\mathrm{training}$, $\\Delta \\ln k_\\mathrm{pdep}$, $\\Delta \\ln k_\\mathrm{family}$, $\\Delta \\ln k_\\mathrm{non-exact}$, and $\\Delta \\ln k_\\mathrm{rule}$.\n",
    "\n",
    "For library, training, and pdep reactions, the kinetic uncertainty is assigned according to their uncertainty type.  For kinetics estimated using RMG's rate rules, the following formula is used to calculate the uncertainty:\n",
    "\n",
    "$$\\Delta \\ln k_\\mathrm{rate\\; rules} = \\Delta\\ln k_\\mathrm{family} + \\log_{10}(N+1) \\left(\\Delta\\ln k_\\mathrm{non-exact}\\right)  + \\sum_{\\mathrm{rule}\\; i} w_i \\Delta \\ln k_{\\mathrm{rule},i}$$\n",
    "\n",
    "where N is the total number of rate rules used and $w_{i}$ is the weight of the rate rule in the averaging scheme for that kinetics estimate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: You must load the database with the same settings which were used to generate the model.\n",
    "#       This includes any thermo or kinetics libraries which were used.\n",
    "uncertainty.load_database(\n",
    "    thermo_libraries=['DFT_QCI_thermo', 'primaryThermoLibrary'],\n",
    "    kinetics_families='default',\n",
    "    reaction_libraries=[],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncertainty.extract_sources_from_model()\n",
    "uncertainty.assign_parameter_uncertainties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first order local uncertainty, or variance $(d\\ln c_i)^2$, for the concentration of species $i$ is defined as:\n",
    "\n",
    "$$(\\Delta \\ln c_i)^2 = \\sum_{\\mathrm{reactions}\\; m} \\left(\\frac{\\partial\\ln c_i}{\\partial\\ln k_m}\\right)^2 (\\Delta \\ln k_m)^2  + \\sum_{\\mathrm{species}\\; n} \\left(\\frac{\\partial\\ln c_i}{\\partial G_n}\\right)^2(\\Delta G_n)^2$$\n",
    "\n",
    "We have previously performed the sensitivity analysis.  Now we perform the local uncertainty analysis and apply the formula above using the parameter uncertainties and plot the results.  This first analysis considers the parameters to be independent.  In other words, even when multiple species thermochemistries depend on a single thermo group or multiple reaction rate coefficients depend on a particular rate rule, each value is considered independent of each other.  This typically results in a much larger uncertainty value than in reality due to cancellation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "result = uncertainty.local_analysis(sensitive_species, correlated=False, number=5, fileformat='.png')\n",
    "print(process_local_results(result, sensitive_species, number=5)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Show the uncertainty plots\n",
    "for species in sensitive_species:\n",
    "    print('{}: Thermo Uncertainty Contributions'.format(species))\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory, 'uncorrelated', 'thermoLocalUncertainty_{}.png'.format(species.to_chemkin()))))\n",
    "    \n",
    "    print('{}: Reaction Uncertainty Contributions'.format(species))\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory, 'uncorrelated', 'kineticsLocalUncertainty_{}.png'.format(species.to_chemkin()))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Uncertainty assignment and propagation of correlated parameters\n",
    "\n",
    "A more accurate picture of the uncertainty in mechanism estimated using groups and rate rules requires accounting of the correlated errors resulting from using the same groups in multiple parameters.  This requires us to track the original sources: the groups and the rate rules, which constitute each parameter.  These errors may cancel in the final uncertainty calculation.  Note, however, that the error stemming from the estimation method itself do not cancel.  \n",
    "\n",
    "For thermochemistry, the error terms described previously are $\\Delta G_\\mathrm{library}$, $\\Delta G_\\mathrm{QM}$, $\\Delta G_\\mathrm{GAV}$, and $\\Delta _\\mathrm{group}$.  Of these, $\\Delta G_\\mathrm{GAV}$ is an uncorrelated residual error, whereas the other terms are correlated. The set of correlated and uncorrelated parameters can be thought of instead as a set of independent parameters, $\\Delta G_{ind,w}$.\n",
    "\n",
    "For kinetics, the error terms described perviously are $\\Delta \\ln k_\\mathrm{library}$, $\\Delta \\ln k_\\mathrm{training}$, $\\Delta \\ln k_\\mathrm{pdep}$, $\\Delta \\ln k_\\mathrm{family}$, $\\Delta \\ln k_\\mathrm{non-exact}$, and $\\Delta \\ln k_\\mathrm{rule}$.  Of these, $\\Delta \\ln k_\\mathrm{family}$ and $\\Delta \\ln k_\\mathrm{non-exact}$ are uncorrelated error terms resulting from the method of estimation.  Again, we consider the set of correlated and uncorrelated parameters as the set of independent parameters, $\\Delta\\ln k_{ind,v}$.\n",
    "\n",
    "The first order local uncertainty, or variance $(d\\ln c_i)^2$, for the concentration of species $i$ becomes:\n",
    "\n",
    "$$(\\Delta \\ln c_i)^2 = \\sum_v \\left(\\frac{\\partial\\ln c_i}{\\partial\\ln k_{ind,v}}\\right)^2 \\left(\\Delta\\ln k_{ind,v}\\right)^2 + \\sum_w \\left(\\frac{\\partial\\ln c_i}{\\partial G_{ind,w}}\\right)^2 \\left(\\Delta G_{ind,w}\\right)^2$$\n",
    "\n",
    "where the differential terms can be computed as:\n",
    "\n",
    "$$\\frac{\\partial\\ln c_i}{\\partial\\ln k_{ind,v}} = \\sum_m \\frac{\\partial\\ln c_i}{\\partial\\ln k_m} \\frac{\\partial\\ln k_m}{\\partial\\ln k_{ind,v}}$$\n",
    "\n",
    "$$\\frac{\\partial\\ln c_i}{\\partial G_{ind,w}} = \\sum_n \\frac{\\partial\\ln c_i}{\\partial G_n} \\frac{\\partial G_n}{\\partial G_{ind,w}}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "uncertainty.assign_parameter_uncertainties(correlated=True)\n",
    "result = uncertainty.local_analysis(sensitive_species, correlated=True, number=10, fileformat='.png')\n",
    "print(process_local_results(result, sensitive_species, number=5)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Show the uncertainty plots\n",
    "for species in sensitive_species:\n",
    "    print('{}: Thermo Uncertainty Contributions'.format(species))\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory, 'correlated', 'thermoLocalUncertainty_{}.png'.format(species.to_chemkin()))))\n",
    "    \n",
    "    print('{}: Reaction Uncertainty Contributions'.format(species))\n",
    "    display(Image(filename=os.path.join(uncertainty.output_directory, 'correlated', 'kineticsLocalUncertainty_{}.png'.format(species.to_chemkin()))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:rmg_env]",
   "language": "python",
   "name": "conda-env-rmg_env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
